Introduction à la statistique sous R
Mise à niveau M2 Géoprisme 2021
INITIATION A R-BASE
Cette initiation est destinée à des étudiants n’ayant jamais utilisé R mais connaissant à peu près la statistique univariée et bivariée…
L’idée pédagogique est d’apprendre directement aux étudiants à programmer en R markdown plutôt qu’en R. Pourquoi ? Parce qu’ainsi ils vont simultanément :
- taper du code R qu’ils ignorent
- écrire sous ce code les explications du point de vue informatique
- observer les résultats statistiques
- interpréter ces résultats d’un point de vue statistique
Cela n’a l’air de rien, mais en procédant ainsi les étudiants apprennent à produire à la fois leurs notes de cours en R, leurs notes de cours en statistiques et … le langage Rmarkdown pour rédiger leurs futurs travaux.
Bref, si tout a bien marché, l’étudiant n’aura même pas besoin de consulter le présent document, si ce n’est pour vérifier que son programme donne les mêmes résultats …
Les deux exercices qui suivent utilisent volontairement les fonctions de base du langage R (on dit que l’on programme en R-base) à l’exclusion de tout package c’est-à-dire de tout outil graphique ou statistique mis au point ultérieurement.
Par comparaison avec le jeu de lego, cela revient à effectuer des constructions avec la boîte de base. A première vue cela peut sembler frustrant. Mais en réalité cela ne bride en rien l’imagination et permet d’apprendre plein de choses sans être distrait …
1 R-BASE : Projet, Programme R, Document Rmd
- Au commencement, Dieu créa le ciel et la terre.
- Or la terre était vide et vague, les ténèbres couvraient l’abîme, un vent de Dieu tournoyait sur les eaux.
- Dieu dit : Que la lumière soit et la lumière fut.
1.1 Projet
Si l’on veut s’épargner bien des désagréments dans l’apprentissage de R, il faut prendre dès le départ de bonnes habitudes. Parmi celles-ci, l’une des plus importantes est le fait d’inscrire toujours son travail dans le cadre d’un projet R c’est-à-dire - en simplifiant beaucoup - un répertoire de travail contenant l’ensemble des données, programmes, résultats… que l’on pourra par la suite compresser, archiver et transmettre à quelqu’un d’autre.
1.1.1 Lancement de R studio
Sauf à être complètement masochiste, on n’utilise jamais R directement mais on lance d’abord l’interface R-Studio qui facilite conisdérablement l’ensemble des opérations et offre une gamme considérable de services. Il ne faut toutefois pas confondre les deux et il serait par exemple ridicule d’indiquer sur un CV en vue d’un emploi de statisticien que l’on sait utiliser R-studio en oubliant de préciser que l’on maîtrise R.
1.1.2 Création d’un projet
Pour créer un projet on utilise le menus déroulant File/new project/ … et on définit un dossier de notre ordinateur (existant ou à créer) qui contiendra le projet. Une fois l’opération effectuée, on pourra constater que ce dossier contient un fichier xxx.Rproj ou xxx est en principe le nom du dossier dans lequel vous avez stocké le projet.
Ce fichier contient toute une série d’informations dont nous ne parlerons pas dans ce cours d’initiation mais qui, pour faire simple, définissent un ensemble de réglages du logiciel et de préférences de l’utilisateur.
Si vous fermez Rstudio (faites-le !) il vous suffira pour reprendre votre travail là où vous vous étiez arrêté :
- de lancer R-studio et de cliquer sur File/open project/… suivi du nom du fichier xxx.Rproj
- ou plus simplement encore de double cliquer sur le fichier xxx.Rproj ce qui lancera automatiquement Rstudio
Le dossier contenant votre projet R peut être organisé à votre convenance. Certains mettent tout les fichier pêle-mêle dans le dossier. D’autres préfèrent créer des sous-dossiers contenant des données, des programmes, des résultats, des figures. Vous déciderez à l’usage de ce qui vous convient le mieux, mais le point important est que tout ce qui entre ou sort de vos programmes R doit être de préférence stocké dans le répertoire du projet.
knitr::include_graphics("figures/monprojet.png")1.2 Programme
La fonction initiale d’un langage de programmation comme R est … de créer des programmes c’est-à-dire des ensembles d’instruction permettant d’accomplir une tâche à l’intérieur d’une chaîne de production. Dans le cas d’un logiciel spécialisé dans l’analyse statistique, il s’agira donc de partir de données (statistiques, géographiques, textuelles, …) pour aboutir à des résultats prenant la forme de tableaux, cartes ou graphiques. Il ne s’agit donc en somme que d’une étape du travail de recherche où le principal avantage de R est d’automatiser une tâche et de faciliter sa reproduction ultérieure avec en arrière plan un objectif de productivité puisque l’ordinateur réalise en quelques millisecondes des tâches qui prendraient des heures avec un logiciel click-bouton de type Excel.
knitr::include_graphics("figures/ProgrammeR.png")Prenons un exemple simple de problème facile à résoudre avec R mais plus compliqué avec des logiciels click-boutons. Il s’agit d’un exemple pédagogique tiré d’un très vieux cours d’analyse spatiale portant sur les semis de point et les localisations optimales.
On considère 5 station services localisés à l’intérieur d’une ville à plan en damier qui livrent chacune la même quantité de carburant par semaine aux clients. Où faut-il localiser les deux infrastructures suivantes : - un dépôt de carburant permettant d’alimenter les cinq stations qui doit minimiser la distance moyenne de livraison. - une caserne de pompier qui doit pouvoir intervenir rapidement sur toute les stations et qui doit minimiser la distance maximale à la station la plus éloignée.
knitr::include_graphics("figures/access.png")On constitue deux équipes d’étudiants, certains utilisant un programme R et d’autres Excel. On se propose de voir qui ira le plus vite :
1.2.1 Round 1. Saisie des données et affichage du tableau
On crée un programme R avec File/New File/R Script puis on l’enregistre avec File/Save/ … suivi du nom du programme.
# Saisie des variables
CODE <- c("A","B","C","D","E")
X <- c(10,20,40,50,180)
Y <- c(40,60,40,60,50)
# Regroupement dans un tableau
coo <- data.frame(X,Y)
# Ajout du nom des lignes
row.names(coo) <- CODE
# Affichage du tableau
coo X Y
A 10 40
B 20 60
C 40 40
D 50 60
E 180 50
Normalement, les étudiants qui utilisent un tableur ont du aller plus vite et Excel mène sur R par 1-0
1.2.2 Round 2. Affichage de la carte
Vous devez essayez de reproduire l’image correspondant au problème posé
plot(X,Y,
col="red",
pch=20,
xlim=c(0,180),
ylim=c(0,90),
asp = 1)
text(X,Y,
labels = CODE,
pos = 2) La création d’un graphique est à première vue plus facile avec un logiciel click-bouton. L’avantage est très clairement pour Excel qui mène désormais 2 à 0.
1.2.3 Round 3. Calcul de la station la plus accessible à vol d’oiseau (distance euclidienne)
Vous devez calculer une matrice de distance euclidienne entre toutes les stations et trouver la plus accessible.
# calcul la matrice de distance euclidienne
mat<-dist(coo, upper = T, method = "euclidean")
mat A B C D E
A 22.36068 30.00000 44.72136 170.29386
B 22.36068 28.28427 30.00000 160.31220
C 30.00000 28.28427 22.36068 140.35669
D 44.72136 30.00000 22.36068 130.38405
E 170.29386 160.31220 140.35669 130.38405
# distance moyenne
apply(as.matrix(mat),1,mean) A B C D E
53.47518 48.19143 44.20033 45.49322 120.26936
Là, je parie que les utilisateurs d’Excel ont eu un peu plus de mal … En tous les cas, Excel ne mèneplus que par 2 à 1
1.2.4 Round 4. Calcul de la station la plus accessible par la route (distance de Manhattan)
Vous devez calculer une matrice de distance de Manhattan entre toutes les stations et trouver la plus accessible.
# calcul la matrice de distance de Manhattan
mat<-dist(coo,upper = T, method = "manhattan")
mat A B C D E
A 30 30 60 180
B 30 40 30 170
C 30 40 30 150
D 60 30 30 140
E 180 170 150 140
# distance moyenne de Manhattan
apply(as.matrix(mat),1,mean) A B C D E
60 54 50 52 128
Je reconnais que c’est unpeu facile, mais à nouveau R l’emporte ce qui fait désormais match nul 2-2
1.2.5 Round 5. Localisation du dépôt de carburant
Dans le cas particulier de la distance de Manhattan, le calcul du point le plus proche de tous les autres s’obtient facilement en calculant le point médian dont les coordonnées correspondent à la médiane de X et la médiane de Y.
medX <- median(X)
medX[1] 40
medY <- median(Y)
medY[1] 50
A priori, le calcul est aussi facile dans R et dans Excel : match nul 3-3
1.2.6 Round 6. Localisation de la caserne de pompiers
Dans le cas particulier de la distance de Manhattan, le calcul du point minimisant la distance maximale s’obtient en trouvant le centre du diamètre minimal en X et en Y. Il s’agit de la localisation la plus équitable où le plus défavorisé est le moins défavorisé possible.
equX <- (max(X)+min(X))/2
equX[1] 95
equY <- (max(Y)+min(Y))/2
equY[1] 50
A priori, le calcul est toujours aussi facile dans R et dans Excel : match nul 4-4
1.2.7 Round 7. Visualisation des deux points sur la carte
On va placer en bleu le point médian et en vert le point le plus équitable. Dans le cas de R on peut recopier les lignes de code du graphique du round n°2 ce qui gagne désormais du temps :
# Programme antérieur
plot(X,Y,
col="red",
pch=20,
xlim=c(0,180),
ylim=c(0,90),
asp = 1)
text(X,Y,
labels = CODE,
pos = 2)
# Ajout du dépôt de carburant
points(medX, medY, col="blue", pch=3)
text(medX,medY, "dépot",pos=1)
# Ajout du point médian
points(equX, equY, col="green", pch=3)
text(equX,equY, "caserne",pos=1)Le résultat du match est incertain mais R n’est plus désavantagé puisqu’on peut recycler les lignes de code précédentes pour le graphique de base. Disons 5-5 même s’il y a de fortes chances que R l’emporte.
1.2.8 Dernier round. Refaire toute l’analyse avec une station de plus
Deux stations F(100,20) et G(150,30) ont été ajoutées. Il faut refaire la carte finale. Cela ne pose aucun problème dans R puisqu’il suffit de modifier l’entrée des données et récupérer des bouts de programme
# (1) Saisie des variables
CODE <- c("A","B","C","D","E","F","G")
X <- c(10,20,40,50,180,100,150)
Y <- c(40,60,40,60,50,20,30)
coo <- data.frame(X,Y)
row.names(coo) <- CODE
# (2) calcul des points centraux
medX <- median(X)
medY <- median(Y)
equX <- (max(X)+min(X))/2
equY <- (max(Y)+min(Y))/2
# (3) Graphique
plot(X,Y,
col="red",
pch=20,
xlim=c(0,180),
ylim=c(0,90),
asp = 1)
text(X,Y,
labels = CODE,
pos = 2)
# Ajout du dépôt de carburant
points(medX, medY, col="blue", pch=3)
text(medX,medY, "Dépôt",pos=1)
# Ajout du point médian
points(equX, equY, col="green", pch=3)
text(equX,equY, "Caserne",pos=1)Excel n’a aucune chance d’aller plus vite et R remporte le match par KO !
1.3 Document R markdown
knitr::include_graphics("figures/DocumentRmd.png")2 R-BASE : Variables quantitatives, corrélation, régression, test d’égalité des moyennes
Accrochez vos ceintures … on plonge directement dans un programme R qu’il s’agit juste d’exécuter sans comprendre ce qui se passe pour avoir une idée des possibilités du logiciel mais aussi des difficultés. Vous êtes un peu comme un étranger qui entend parler une langue nouvelle et découvre de nouveaux mots et essaye de les reproduire comme dans la méthode Assimil … Evidemment ça ne marche pas à tous les coups.
On commence ici par un exemple où tout se passe bien et on utilise pour cela un exemple fétiche de l’enseignant, celui des pays d’Europe en 1988 à la veille de la chute du mur de Berlin. Exemple qui a été testé depuis 30 ans sur tous les logiciels de statistique possible (Lotus, Calc, SAS, XLSTAT, SPAD, ….). Il vous est recommandé d’avoir également un exemple fétiche sur lequel vous pourrez tester à votre tour le programme ci-dessous en l’adaptant.
Les 24 pays d’Europe sont identifiés par leur code iso (PAYS), leur appartenance aux pays socialistes ou capitalistes (BLOC), leur position spatiale (X,Y) et tout un ensemble de variables économpiques et sociale telles que le PNB par habitant (PNB), le taux de mortlaité infantile (TMI), l’espérance de vie à la naissance (ESP), le taux d’urbanisation (URB), le taux de natalité (NAT), le taux de mortalité (MOR), l’indice conjoncturel de fécondité (FEC), le % de jeunes de 0-14 ans (JEU), le % de vieux de 65 ans et plus (VIE), la superficie totale en milliers de kM2 (SUP) et la population totale en millions (POP).
2.1 Définir un environnement de travail
2.1.1 Définir un espace de travail
le répertoire de travail dans lequel seront lus ou écrit les fichiers.
Dans le programme ci-dessous, les lignes commençant par “#” sont des commentaires. La seule ligne active est getwd() qui indique la position du répertoire courant. Il est recommandé aux débutants, mais aussi aux programmeurs plus expérimentés de ne pas hésiter à ajouter beaucoup de lignes de commentaires. Cela facilité l’apprentissage et permet de transmettre ses programmes à d’autres.
<<<<<<< HEAD# ---------------------------
# (1) ESPACE DE TRAVAIL
# ---------------------------
# (1.1) Quel est le répertoire actuel ?
getwd()[1] "/Users/claudegrasland1/git/Init_R_2021"
=======
# ---------------------------
# (1) ESPACE DE TRAVAIL
# ---------------------------
# (1.1) Quel est le répertoire actuel ?
getwd()[1] "C:/git/Init_R_2021"
>>>>>>> d676e259c8be1e9f88c818dd0e1c0cc0a591abdd
On repère ensuite le chemin de l’emplacement du dossier où se trouvent les données et on examine son contenu. Attention, contrairement à windows il faut toujours utiliser des slash (/) et non pas des antislash pour décrire les chemins d’accès des fichiers
# (1.2) Choix du rÈpertoire contenant les donnÈes
monchemin<-"data/euro1988"
list.files(monchemin)[1] "euro1988.rda" "euro1988.txt"
Nous allons maintenant essayer de charger le fichier de données
2.2 Préparer son tableau de données
On utilise ici une procédure adaptée à la lecture de fichiers texte : read.table() L’instruction header=TRUE signale que la première ligne donne le nom des variables L’instruction dec=“.” signale que les décimales sont représentées par des points et non pas des virgules.
2.2.1 Importation d’un tableau
# --------------------------------------------------------
# (2) IMPORTATION ET MISE EN FORME D'UN TABLEAU DE DONNEES
# --------------------------------------------------------
# (2.1) Importation d'un fichier .txt
euro <- read.table("data/euro1988/euro1988.txt",
dec=".",
header=TRUE)Pour savoir si le chargement s’est bien pssé on peut afficher le tableau en tapant son nom
euro PAYS BLOC PNB TMI ESP URB NAT MOR FEC JEU VIE SUP POP X Y
1 ALB Soc 600 43.0 71 34 27 6 3.3 35 5 29 3.1 1445 684
2 AUT Cap 10000 10.3 75 55 12 12 1.4 18 14 84 7.6 1121 1031
3 BEL Cap 9200 9.7 75 95 12 11 1.5 19 14 31 9.9 734 1187
4 BUL Soc 2000 14.5 72 65 13 11 2.0 21 11 111 9.0 1649 839
5 DAN Cap 12600 8.4 75 84 11 11 1.5 18 15 43 5.1 887 1500
[ reached 'max' / getOption("max.print") -- omitted 19 rows ]
Mais on peut aussi se contenter d’afficher les premières lignes (head) ou les dernières lignes (tail) du tableau. Par exemple, pour voir les trois premières et les 5 dernières lignes :
head(euro,3) PAYS BLOC PNB TMI ESP URB NAT MOR FEC JEU VIE SUP POP X Y
1 ALB Soc 600 43.0 71 34 27 6 3.3 35 5 29 3.1 1445 684
2 AUT Cap 10000 10.3 75 55 12 12 1.4 18 14 84 7.6 1121 1031
3 BEL Cap 9200 9.7 75 95 12 11 1.5 19 14 31 9.9 734 1187
tail(euro,5) PAYS BLOC PNB TMI ESP URB NAT MOR FEC JEU VIE SUP POP X Y
20 ROY Cap 8900 9.5 75 91 13 12 1.8 19 15 245 57.1 336 1452
21 SUE Cap 13200 5.9 77 83 12 11 1.8 18 18 450 8.4 1061 1931
22 SUI Cap 17800 6.8 77 61 12 9 1.5 17 14 41 6.6 882 956
23 TCH Soc 3200 13.9 71 74 14 12 2.0 24 11 128 15.6 1219 1157
24 YOU Soc 2300 27.1 70 47 15 8 2.1 24 8 256 23.6 1356 855
Naturellement … les instructions présentées ci-dessus ne marcheront pas pour tous les tableaux et vous devez vous attendre à pas mal de difficultés à ce stade. Il existe en effet beaucoup d’options à connaître pour charger les fichiers. Vous pouvez ainsi examiner toutes les options de la procédure read.table en tapant un ? suivi du nom de la procédure. Et vous pouvez faire de même pour les autres procédures d’importation comme read.csv, read.delim, etc…
2.2.2 Choix du nom des lignes
Par défaut, le nom des lignes correspond à l’ordre de lecture du fichier (1..24). Mais il peut être intéressant de le préciser en lui attribuant par exemple la valeur d’une variable servant d’identifiant.
rownames(euro)<-euro$PAYS
head(euro) PAYS BLOC PNB TMI ESP URB NAT MOR FEC JEU VIE SUP POP X Y
ALB ALB Soc 600 43.0 71 34 27 6 3.3 35 5 29 3.1 1445 684
AUT AUT Cap 10000 10.3 75 55 12 12 1.4 18 14 84 7.6 1121 1031
BEL BEL Cap 9200 9.7 75 95 12 11 1.5 19 14 31 9.9 734 1187
BUL BUL Soc 2000 14.5 72 65 13 11 2.0 21 11 111 9.0 1649 839
DAN DAN Cap 12600 8.4 75 84 11 11 1.5 18 15 43 5.1 887 1500
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
2.2.3 Vérification du type des variables
On peut s’éviter beaucoup d’ennuis en contrôlant le type des variables qui a été attribué aux différentes colonnes par R lors de la lecture d’un fichier. Pour cela on commence par regarder quels sont les types de variables du tableau à l’aide de l’instruction str() :
str(euro)'data.frame': 24 obs. of 15 variables:
$ PAYS: chr "ALB" "AUT" "BEL" "BUL" ...
$ BLOC: chr "Soc" "Cap" "Cap" "Soc" ...
$ PNB : int 600 10000 9200 2000 12600 4800 12200 10100 3700 2000 ...
$ TMI : num 43 10.3 9.7 14.5 8.4 9 5.8 8 12.3 19 ...
$ ESP : int 71 75 75 72 75 76 74 75 74 70 ...
$ URB : int 34 55 95 65 84 91 62 73 58 58 ...
$ NAT : int 27 12 12 13 11 12 12 14 11 12 ...
$ MOR : int 6 12 11 11 11 8 10 10 9 14 ...
$ FEC : num 3.3 1.4 1.5 2 1.5 1.7 1.6 1.8 1.7 1.8 ...
$ JEU : int 35 18 19 21 18 23 19 21 21 21 ...
$ VIE : int 5 14 14 11 15 12 13 13 13 13 ...
$ SUP : int 29 84 31 111 43 505 337 551 132 93 ...
$ POP : num 3.1 7.6 9.9 9 5.1 39 4.9 55.9 10.1 10.6 ...
$ X : int 1445 1121 734 1649 887 358 1299 646 1453 1341 ...
$ Y : int 684 1031 1187 839 1500 582 2103 941 593 1042 ...
L’instruction str() est très importante dans R. Elle permet en effet d’examiner le type d’un objet et des éléments qui le composent. Ici on note que l’objet euro est un data.frame c’est-à-dire un tableau de données croisant des indivdus (lignes) et des variables de types différents (colonnes). Par défaut, R a attribué le type “Factor” à toutes les variables composées de texte, le type “int” ou integer à toutes les variables comportant uniquement des nombres entiers et enfin le type “num” ou numérique à toutes les variables comportant des chiffres avec des décimales. Il peut arriver que ce choix par défaut ne correspondent pas aux traitements que l’on veut effectuer et, dans ce cas, on procèdera à des conversions de variables d’un type à un autre.
Ainsi le type “Factor” est difficile à utiliser pour les débutants et il est souvent confondu avec le type “Character”. Pour fixer les idées, on peut dire qu’un “Factor” correspond à ces classes portant à la fois un numéro et une étiquette. Ainsi la variable “BLOC” est bien un facteur correspondant à deux classes : celle des pays capitalistes (numéro = 1, étiquette =“cap”) et celle des pays socialistes (numero = 2, étiquette = “soc”). Par contre la variable PAYS ne définit pas des classes puisqu’elle est différente pour chaque pays et est donc du type “character”.
De la même manière, un nombre entier peut conduire à des comportements bizarre pour les valeurs élevées et il vaut mieux utiliser le tupe numérique lorsqu’on est en présence de nombre réels qui ont juste été arrondis à zéro chiffres après la virgule. Dans notre exemple, il n’existe pas de variables entières mais uniquement des nombres réels arrondis.
euro$PAYS<-as.character(euro$PAYS)
euro$BLOC<-as.factor(euro$BLOC)
euro$PNB<-as.numeric(euro$PNB)
euro$ESP<-as.numeric(euro$ESP)
euro$URB<-as.numeric(euro$URB)
euro$NAT<-as.numeric(euro$NAT)
euro$MOR<-as.numeric(euro$MOR)
euro$JEU<-as.numeric(euro$JEU)
euro$VIE<-as.numeric(euro$VIE)
euro$SUP<-as.numeric(euro$SUP)
euro$POP<-as.numeric(euro$POP)
euro$X<-as.numeric(euro$X)
euro$Y<-as.numeric(euro$Y)
str(euro)'data.frame': 24 obs. of 15 variables:
$ PAYS: chr "ALB" "AUT" "BEL" "BUL" ...
$ BLOC: Factor w/ 2 levels "Cap","Soc": 2 1 1 2 1 1 1 1 1 2 ...
$ PNB : num 600 10000 9200 2000 12600 4800 12200 10100 3700 2000 ...
$ TMI : num 43 10.3 9.7 14.5 8.4 9 5.8 8 12.3 19 ...
$ ESP : num 71 75 75 72 75 76 74 75 74 70 ...
$ URB : num 34 55 95 65 84 91 62 73 58 58 ...
$ NAT : num 27 12 12 13 11 12 12 14 11 12 ...
$ MOR : num 6 12 11 11 11 8 10 10 9 14 ...
$ FEC : num 3.3 1.4 1.5 2 1.5 1.7 1.6 1.8 1.7 1.8 ...
$ JEU : num 35 18 19 21 18 23 19 21 21 21 ...
$ VIE : num 5 14 14 11 15 12 13 13 13 13 ...
$ SUP : num 29 84 31 111 43 505 337 551 132 93 ...
$ POP : num 3.1 7.6 9.9 9 5.1 39 4.9 55.9 10.1 10.6 ...
$ X : num 1445 1121 734 1649 887 ...
$ Y : num 684 1031 1187 839 1500 ...
Et voilà : désormais notre tableau ne comporte plus que des variables de type caractère ou de type numérique. Cela nous a demandé un peu de temps mais nous évitera par la suite beaucoup d’ennuis. En effet, au moins 50 % des erreurs dans R viennent d’un problème de mauvaise définiton du type d’’un objet ou du type d’une variable !!!
2.2.4 Enregistrement du tableau au format propre à R
Il n’est pas forcément necessaire d’enregistrer ses données au format propre à R (.rda)s’il s’agit d’un petit tableau. Il suffit en effet de sauvegarder le programme que nous venons d’écrire et de l’executer chaque fois que nécessaire. Toutefois, il peut arriver que cette étape soit longue et dans ce cas la sauvegarde est utile.
save(euro, file="data/euro1988/euro1988.rda")on pourra ensuite facilement récupérer les données par la commande suivante
load("data/euro1988/euro1988.rda")2.2.5 Un petit résumé statistique
N’oublions pas que R est avant tout un logiciel de statistiques… La récompense d’une bonne définition du tableau sera d’obtenir pour chaque variable un résumé adapté à son type :
summary(euro) PAYS BLOC PNB TMI ESP
Length:24 Cap:16 Min. : 600 Min. : 5.800 Min. :70.00
Class :character Soc: 8 1st Qu.: 2275 1st Qu.: 8.475 1st Qu.:71.75
Mode :character Median : 6850 Median : 9.600 Median :74.50
Mean : 7208 Mean :13.121 Mean :73.67
3rd Qu.:10575 3rd Qu.:14.825 3rd Qu.:75.00
URB NAT MOR FEC
Min. :30.00 Min. :10.00 Min. : 6.00 Min. :1.400
1st Qu.:57.50 1st Qu.:12.00 1st Qu.: 9.00 1st Qu.:1.575
Median :68.00 Median :12.50 Median :10.50 Median :1.700
Mean :67.92 Mean :13.46 Mean :10.33 Mean :1.829
3rd Qu.:83.25 3rd Qu.:14.00 3rd Qu.:11.00 3rd Qu.:2.000
JEU VIE SUP POP
Min. :15.00 Min. : 5.0 Min. : 29.0 Min. : 3.10
1st Qu.:19.00 1st Qu.:11.0 1st Qu.: 80.5 1st Qu.: 7.35
Median :20.00 Median :13.0 Median :130.0 Median :10.45
Mean :21.33 Mean :12.5 Mean :198.5 Mean :20.64
3rd Qu.:23.25 3rd Qu.:14.0 3rd Qu.:304.0 3rd Qu.:27.20
X Y
Min. : 147.0 Min. : 567
1st Qu.: 760.2 1st Qu.: 851
Median :1042.0 Median :1100
Mean : 991.7 Mean :1154
3rd Qu.:1309.5 3rd Qu.:1348
[ getOption("max.print") est atteint -- ligne 1 omise ]
On va tenter de répondre à des questions simples sans détailler dans l’immédiat les outils statistiques qui permettent de valider scientifiquement les réponses. On insistera plus ici sur les capacités de visualisation de résultats sous R.
2.3 Expliquer une variable quantitative (Y) par une variable qualitative (X)
La mortalité infantile (Y) est-elle significativement plus forte dans les pays socialistes (X)?
Notre objectif est plus généralement de mettre en relation une variable numérique (TMI) avec une variable de type facteur (BLOC). On va procéder par étape en étudiant d’abord chque variable séparément puis en les croisant et finalement en testant notre hypothèse d’une surmortalité infantile des pays socialistes.
2.3.1 Analyse de la répartition des pays par bloc politique (BLOC)
Deux petites commandes pour dénombrer l’effectif des classes et les représenter
# Analyse de la variable qualitative (BLOC)
#calcul des fréquences
summary(euro$BLOC)Cap Soc
16 8
<<<<<<< HEAD
#diagramme en bâtons
plot(euro$BLOC)#diagramme en bâtons
plot(euro$BLOC)2.3.2 Analyse de la distribution des taux de mortalité infantile (TMI)
Quelques petites commandes pour camlculer les paramètres principaux de la mortalité infantile et la représenter.
# Analyse de la variable quantitative (TMI)
# calcul et visualisation des quantiles
summary(euro$TMI) Min. 1st Qu. Median Mean 3rd Qu. Max.
5.800 8.475 9.600 13.121 14.825 43.000
<<<<<<< HEAD
boxplot(euro$TMI)#calcul de la moyenne et de l'écart type
mean(euro$TMI)boxplot(euro$TMI)#calcul de la moyenne et de l'écart type
mean(euro$TMI)[1] 13.12083
sd(euro$TMI)[1] 8.513875
<<<<<<< HEAD
# visualisation de l'histogramme
hist(euro$TMI)# visualisation de l'histogramme
hist(euro$TMI)2.3.3 Comaparaison de la mortalité infantile par bloc (BLOC*TMI)
Quelques outils pour analyser conjointement les deux distributions.
# calcul et visualisation des quantiles
tapply(euro$TMI, euro$BLOC, summary)$Cap
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.800 7.925 8.650 9.069 9.800 15.800
$Soc
Min. 1st Qu. Median Mean 3rd Qu. Max.
9.20 14.35 18.25 21.23 25.98 43.00
<<<<<<< HEAD
boxplot(euro$TMI~euro$BLOC)#calcul de la moyenne et de l'écart type
tapply(euro$TMI, euro$BLOC, mean)boxplot(euro$TMI~euro$BLOC)#calcul de la moyenne et de l'écart type
tapply(euro$TMI, euro$BLOC, mean) Cap Soc
9.06875 21.22500
tapply(euro$TMI, euro$BLOC, sd) Cap Soc
2.434945 10.624197
<<<<<<< HEAD
# visualisation des deux histogrammes
par(mfrow= c(2,1))
hist(euro$TMI[euro$BLOC =="Soc"],
# xlim = c(0,50),
main="Pays Socialistes",
col="red")
hist(euro$TMI[euro$BLOC =="Cap"],
# xlim = c(0,50),
main="Pays Capitalistes",
col="blue") # visualisation des deux histogrammes
par(mfrow= c(2,1))
hist(euro$TMI[euro$BLOC =="Soc"],
# xlim = c(0,50),
main="Pays Socialistes",
col="red")
hist(euro$TMI[euro$BLOC =="Cap"],
# xlim = c(0,50),
main="Pays Capitalistes",
col="blue") 2.3.4 Test d’égalité des moyennes (BLOC*TMI)
Maintenant que nous avons pris connaissance de la forme de nos deux variables BLOC et TMI, essayons de réaliser une opération statistique de plus haut niveau : le test d’une hypothèse. Nous savons grâce aux analyses précédentes que la mortalité infantile semble plus élévée en général dans les pays socialistes que dans les pays capitalistes. Mais pouvons nous le *prouver" de façon rigoureuse en utilisant les règles de l’art ?
Bien évidemment, ce n’est pas l’apprentissage de R qui va vous apprendre la théorie des tests statistiques. Mais R va vous permettre très rapidement de mettre en oeuvre les tests que vous aurez appris en cours ou dans les manuels. Et surtout, il va vous offrir une gamme de possibilités de tests et de variantes qui devrait vous pousser à renforcer simultanément vos connaissances en R et en statitiques.
Dans l’exemple qui nous intéresse, les étudiants ayant une formation de base en statistique savent qu’il existe un test d’égalité des moyennes de deux échantillons appelé test t de student qui semble approprié au problème de comapraison de la mortalité infantile des pays socialistes et capitalistes. Il tient compte de l’effectif de ces échantillons (combiens de pays de chaque type ?) et de leur hétérogénéité interne (appelée variance). On devine intuitivement que la différence des moyennes de deux échantillons sera d’autant plus significative qu’elle oppose deux groupes homogènes comportant chacun de nombreux pays.
Si cela ne vous parait pas intuitif, imaginez un premier groupe de trois pays ayant pour mortalité (10, 20, 30) et un second groupe de trois pays également ayant pour mortalité (5,15,25). Leurs moyennes respectives sont 20 et 15, soit un écart de 5. Mais il semblerait hasardeux de conclure à de réelles différences compte tenu de l’étalement de chaque distribution.
<<<<<<< HEAD# Comparaison de la distribution de deux échantillons
VARQUANT<-c(10,20,30,5,15,25)
VARQUALI<-c("A","A","A","B","B","B")
par(mfrow=c(1,1))
boxplot(VARQUANT~ VARQUALI)Imaginez maintenant un premier groupe de 6 pays ayant pour mortalité (18,19,20,20,21,22) et un second groupe de 10 pays ayant pour mortalité (13,14,15,15,16,17). Les moyennes des deux groupes sont toujours respectivement de 20 et de 15, mais on “sent” que les écarts sont plus pertinents car ils portent sur deux groupes plus larges et surtout plus homogènes. Dit autrement, il est peu probable que les différences observées dans la seconde expérience soient le fruit du hasard alors que, dans le premier cas, on pouvait imaginer que les pays des deux groupes n’étaient pas fondamentalement différents.
# Comparaison de la distribution de deux échantillons
VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
par(mfrow=c(1,1))
boxplot(VARQUANT~ VARQUALI)# Comparaison de la distribution de deux échantillons
VARQUANT<-c(10,20,30,5,15,25)
VARQUALI<-c("A","A","A","B","B","B")
par(mfrow=c(1,1))
boxplot(VARQUANT~ VARQUALI)Imaginez maintenant un premier groupe de 6 pays ayant pour mortalité (18,19,20,20,21,22) et un second groupe de 10 pays ayant pour mortalité (13,14,15,15,16,17). Les moyennes des deux groupes sont toujours respectivement de 20 et de 15, mais on “sent” que les écarts sont plus pertinents car ils portent sur deux groupes plus larges et surtout plus homogènes. Dit autrement, il est peu probable que les différences observées dans la seconde expérience soient le fruit du hasard alors que, dans le premier cas, on pouvait imaginer que les pays des deux groupes n’étaient pas fondamentalement différents.
# Comparaison de la distribution de deux échantillons
VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
par(mfrow=c(1,1))
boxplot(VARQUANT~ VARQUALI)Le test T de student dont vous trouverez la formule dans tous les bons manuels va permettre de transformer l’intuition précédente en mesure objective et en donner une mesure fondamentale appelée p-value ou seuil de rejet de l’hypothèse d’indépendance entre deux caractères. Là encore, on peut donner une idée intuitive de ce qu’est la p-value en la présentant comme une mesure de la significativité de la relation entre deux caractères. De façon plus rigoureuse, la p-value mesure la probabilité que la relation observée entre deux variables (ici, la différence de moyenne) soit l’effet du hasard. La p-value varie entre 0 et 1, 0 signifiant que la significativité de la relation est certaine (la différence de moyenne ne peut pas être l’effet du hasard) et 1 signifiant qu’il n’existe aucune relation entre les deux variobles (les deux moyennes sont égales).
Examinons brièvement les valeurs de p-value sur nos deux exemples :
# Test paramétrique d'égalité des moyennes de Student
VARQUANT<-c(10,20,30,5,15,25)
VARQUALI<-c("A","A","A","B","B","B")
t.test(VARQUANT ~ VARQUALI)
Welch Two Sample t-test
data: VARQUANT by VARQUALI
t = 0.61237, df = 4, p-value = 0.5734
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.66958 27.66958
sample estimates:
mean in group A mean in group B
20 15
Vous pouvez lire dans le résultat que la p-value est de 0.5734 ce qui signifie que l’on aurait 57% de chances de se tromper si l’on affirmait que la différence de moyenne entre les deux échantillons A et B est l’effet du hasard. Voyons ce qu’il en est pour notre second exemple
# Test paramétrique d'égalité des moyennes de Student
VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
t.test(VARQUANT ~ VARQUALI)
Welch Two Sample t-test
data: VARQUANT by VARQUALI
t = 6.1237, df = 10, p-value = 0.0001121
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
3.180732 6.819268
sample estimates:
mean in group A mean in group B
20 15
Cette fois-ci, vous pouvez voir que la valuer du t de student est beaucoup plus forte (6.13) ce qui se traduit par une p-value très proche de zéro (0.00011). La probabilité que les différences entre les deux échantillons soit l’effet du hasard est donc minime et on peut conclure avec une quasi-certiude que A et B correpondent à des groupes de pays ayant des mortalité significativement différentes.
Essayons maintenant d’appliquer ce même test de Student à notre exemple empirique de la mortalité infantile des pays socialistes et capitalistes.
# Test paramétrique d'égalité des moyennes de Student (inadapté)
t.test(euro$TMI ~ euro$BLOC)
Welch Two Sample t-test
data: euro$TMI by euro$BLOC
t = -3.1946, df = 7.3701, p-value = 0.01417
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-21.06338 -3.24912
sample estimates:
mean in group Cap mean in group Soc
9.06875 21.22500
Cette fois-ci le t de Student est égal à -3.19 et la p-value est égale à 0.0147, ce qui signifie que nous avons environ 1.5% de chances que les différences observées entre pays socialistes et capitalistes soient le simple effet du hasard et ne reflètent pas une véritable opposition. Du coup, que faut-il conclure ?
En sciences sociales, il est de tradition de considérer comme significatives les relations comportant un risque d’erreur inférieur à 5% soit 0.05. Puisque notre p-value est inférieure à 0.05, nous pouvons donc nous appuyer sur cette tradition et écrire que “l’on peut affirmer avec un risque d’erreur inférieur à 5% que les pays socialistes et capitalistes présentaient des différences significatives de mortalité infantile en 1988”.
Ou bien écrire plus simplement “En 1988 en Europe, la mortalité infantile moyenne des pays socialistes (21.2 P.1000) était significativement plus forte que celles des pays capitalistes (9.1 p.1000)” et ajouter dans une note en bas de page la remarque qui fait chic p-value=0.0147.
Ne concluez cependant pas trop vite que vous maitrisez désormais toutes les subtilités du test statistique d’égalité des moyennes. Si notre raisonnement précédent est approximativement juste, il serait sans doute critiqué par un vrai statisticien qui noterait immédiatement que le test t de student n’est pas le plus adapté à notre problème car (1) les distributions de mortalité ne sont pas gaussiennes et comportent des valeurs exceptionnelles comme l’Albanie pour les pays socialites ou le Portugal pour les pays capitalistes et (2) la variance des deux échantillons (leur hétérogénéité interne) n’est pas égale ce qui peut fausser l’usage du test de Student. Dans le cas présent, un bon statisticien préférerait sans doute utiliser un test non paramétrique, fondé sur la comparaison de l’ordre des observations et non pas leur valeur.
Soit par exemple le test de Wilcoxon, qu’il est très facile de mettre en oeuvre sour R en changeant à peine notre programme.
# Test non paramétrique d'égalité des médianes de Wilcoxon (correct)
wilcox.test(euro$TMI~euro$BLOC)
Wilcoxon rank sum exact test
data: euro$TMI by euro$BLOC
W = 8, p-value = 0.0001822
alternative hypothesis: true location shift is not equal to 0
LA p-value est désormais égale à 0.0001 et on peut conclure avec certitude à l’existence de différence entre nos deux échantillons, ce qui n’était pas le cas avec le test de student, moins adapté au cas étudié.
2.4 Expliquer une variable quantitative (Y) par une variable quantitative (X)
La mortalité infantile (Y) peut-elle être estimée à l’aide du PIB/habitant (X) ?
Cette exemple de mise en relation de deux variables quantitatives permet de passer rapidement en revue les principe de base de l’analyse des corrélations entre plusieurs variables et de la réalisation d’un modèle d’ajustement linéaire. Il faut toutefois noter que l’exemple a surtout des vertus pédagogiques car le PNB/hab. des pays socialistes d’Europe en 1988 n’était pas véritablement connu et reposait sur des estimations.
2.4.1 Calcul des coefficients de corrélations
La façon la plus simple de calculer un coefficient de corrélation est d’appliquerla fonction cor à un tableau ne contenant que des variables quantitatives
# Coefficients de corrélation de pearson
tab<-euro[,c("TMI","PNB","URB","FEC","POP")]
cor(tab) TMI PNB URB FEC POP
TMI 1.0000000 -0.67645013 -0.6523693 0.8141417 -0.13586790
PNB -0.6764501 1.00000000 0.4929530 -0.6045034 0.01784716
URB -0.6523693 0.49295300 1.0000000 -0.5436824 0.39790430
FEC 0.8141417 -0.60450337 -0.5436824 1.0000000 -0.20123367
POP -0.1358679 0.01784716 0.3979043 -0.2012337 1.00000000
R calcule par débaut le coefficient de corrélation linéaire de Bravais-Pearson qui varie entre -1 (corrélation négative parfaite) et +1 (corrélation positive parfaite). On peut ainsi voir que le taux de mortalité infantile (TMI) a une très forte corrélation positive avec l’indice conjonctuel de fécondité (+0.81) et une forte corrélation négative avec le PNB/hab (-0.68) ou le taux d’urbanisation (-0.65). La corrélationavec la population du pays est négative mais faible (-0.14).
2.4.2 Test de significativité d’un coefficient de corrélation
Cette procédure n’est cependant pas totalement convaincante d’un point de vue statistique car elle n’affiche pas le degré de significativité de la liaison statistique, la fameuse p-value dont nous avons discuté dans l’exemple précédent. Cette significativité dépend en effet à la fois (1) de la valeur absolue du coefficient de corrélation (plus on se rapproche de +1 ou -1,plus la relation est significative) mais aussi (2) de la taille de l’échantillon utilisé pour calculer cette corrélation, appelée nombre de degrés de liberté. Lorsqu’on calcule un coefficient de corrélation, le nombre de dégré de liberté est égal à (N-2) c’est-àdire au nombre d’observations (N) moins le nombre de paramètres utilisés pour faire le calcul (soit 2, puisqu’on a besoin de la moyenne de X et de celle de Y).
Pour savoir si une relation est significative, il faut procéder à un test statistique à l’aide de la fonction cor.test. On va par exemple tester la significativité des relations entre TMI et les quatre autres variables.
# Test du coefficient de corrélation de Peason
cor.test(tab$TMI,tab$PNB, type = "pearson")
Pearson's product-moment correlation
data: tab$TMI and tab$PNB
t = -4.3081, df = 22, p-value = 0.0002843
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.8483508 -0.3755262
sample estimates:
cor
-0.6764501
cor.test(tab$TMI,tab$URB, type = "pearson")
Pearson's product-moment correlation
data: tab$TMI and tab$URB
t = -4.0373, df = 22, p-value = 0.0005507
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.835811 -0.337894
sample estimates:
cor
-0.6523693
cor.test(tab$TMI,tab$FEC, type = "pearson")
Pearson's product-moment correlation
data: tab$TMI and tab$FEC
t = 6.5763, df = 22, p-value = 1.296e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6116118 0.9165298
sample estimates:
cor
0.8141417
cor.test(tab$TMI,tab$POP, type = "pearson")
Pearson's product-moment correlation
data: tab$TMI and tab$POP
t = -0.64324, df = 22, p-value = 0.5267
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.511244 0.283042
sample estimates:
cor
-0.1358679
On voit ainsi apparaître deux tableaux, l’un avec les coefficients de corrélation des variables et l’autre avec les p-value associées à chacun de ces coefficients. Il apparaît alors que les corrélations du taux de mortalité infantile sont très significatives (<0.001) avec le PNB/hab, le taux d’urbanisation et l’indice conjoncturel de fécondité. Mais que la corrélation entre mortalité infantile et population est non-significative.
Ajoutons maintenant une autre précaution statistique consistant à calculer un autre coefficient de corrélation, le coefficient de corrélation de rang aussi appelé coefficient de Spearman. L’intérêt de ce coefficient est tout d’abord de limiter le jeu des valeurs exceptionelles (outlier) qui sont liées à la position excentré d’une seule observation. Mais c’est aussi plus généralement de mettre en évidence des relations croissantes ou décroissantes de forme non-linéaire entre deux variables. En présence de ces deux types de singularité, le coefficient de Spearman affiche souvent des valeurs très divergentes de celui de Pearson. Et la découverte d’une discordance entre les deux coefficients est un indice à ne pas négliger.
# Coefficient de corrélation de Spearman
cor(tab, method = "spearman") TMI PNB URB FEC POP
TMI 1.0000000 -0.9018927 -0.5690669 0.50791428 0.16782609
PNB -0.9018927 1.0000000 0.5019591 -0.64715761 -0.18533830
URB -0.5690669 0.5019591 1.0000000 -0.43758300 0.35196875
FEC 0.5079143 -0.6471576 -0.4375830 1.00000000 -0.07105537
POP 0.1678261 -0.1853383 0.3519688 -0.07105537 1.00000000
# Test du cofficient de Spearman
cor.test(tab$TMI,tab$PNB, method = "spearman")
Spearman's rank correlation rho
data: tab$TMI and tab$PNB
S = 4374.4, p-value = 1.763e-09
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.9018927
cor.test(tab$TMI,tab$URB, method = "spearman")
Spearman's rank correlation rho
data: tab$TMI and tab$URB
S = 3608.9, p-value = 0.003707
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.5690669
cor.test(tab$TMI,tab$FEC, method = "spearman")
Spearman's rank correlation rho
data: tab$TMI and tab$FEC
S = 1131.8, p-value = 0.01128
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.5079143
cor.test(tab$TMI,tab$POP, method = "spearman")
Spearman's rank correlation rho
data: tab$TMI and tab$POP
S = 1914, p-value = 0.4314
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.1678261
Dans l’exemple (ancien mais pédagogique) qui a été choisi, il y a clairement une singularité dans la relation entre TMI et PNB puisque la corrélation qui était de (+0.68) avec le coefficient de Pearson grimpe à (+0.90) avec celui de Spearman. Dans le même temps, la relation entre TMI et FEC chute de (+0.81) à (+0.51) ce qui indique une autre singularité. Il faudra donc examiner soigneusement les nuages de points pour comprendre d’où viennent ces divergences, ce que l’on va faire au cours de l’étape suivante.
2.4.3 Visualisation de la forme de la relation entre deux variables quantitatives
On peut visualiser une relations entre deux variables X et Y à l’aide de l’instruction générique plot(). Soit par exemple la relation entre urbanisation et mortalité infantile :
<<<<<<< HEADplot(euro$URB,euro$TMI)On peut également ranger plusieurs graphiques en une seule figure en utilisant l’instruction générale par qui définit le format des graphiques, associée au paramètre mfrow qui indique combien de lignes et de colonnes on souhaite générer. Chacune des case de ce tableau servira successivement à afficher les graphiques demandés par l’instruction plot. Ici, nous déclarons que nous voulons ordonner les graphiques sur 2 lignes et 2 colonnes puis nous executons 4 fois de suite l’instruction plot.
# Un seul graphique de synthèse
par(mfrow= c(2,2))
plot(euro$PNB,euro$TMI)
plot(euro$URB,euro$TMI)
plot(euro$FEC,euro$TMI)
plot(euro$POP,euro$TMI)Ces figures permettent (avec un peu d’habitude …) d’expliquer les divergences observées entre les coefficients de corrélation de Pearson et de Spearman.
par(mfrow= c(1,2))
X<-euro$PNB
Y<-euro$TMI
corlin<-round(cor(X,Y),3)
reglin<-lm(Y~X)
plot(X,Y, main=paste("Coeff Pearson = ",corlin))
abline(reglin,col="red")
rankX<-rank(X)
rankY<-rank(Y)
corrank<-round(cor(rankX,rankY),3)
regrank<-lm(rankY~rankX)
plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
abline(regrank,col="red")- la relation entre TMI et PNB est clairement non linéaire et forme un arc de parabole plutôt qu’une droite. C’est la raison pour laquelle le coefficient de Spearman était beaucoup plus élevé que celui de Pearson. En effet, lorsque l’on transforme les variables en rang, la courbure disparaît.
par(mfrow= c(1,2))
X<-euro$FEC
Y<-euro$TMI
corlin<-round(cor(X,Y),3)
reglin<-lm(Y~X)
plot(X,Y, main=paste("Coeff Pearson = ",corlin))
abline(reglin,col="red")
rankX<-rank(X)
rankY<-rank(Y)
corrank<-round(cor(rankX,rankY),3)
regrank<-lm(rankY~rankX)
plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
abline(regrank,col="red")plot(euro$URB,euro$TMI)On peut également ranger plusieurs graphiques en une seule figure en utilisant l’instruction générale par qui définit le format des graphiques, associée au paramètre mfrow qui indique combien de lignes et de colonnes on souhaite générer. Chacune des case de ce tableau servira successivement à afficher les graphiques demandés par l’instruction plot. Ici, nous déclarons que nous voulons ordonner les graphiques sur 2 lignes et 2 colonnes puis nous executons 4 fois de suite l’instruction plot.
# Un seul graphique de synthèse
par(mfrow= c(2,2))
plot(euro$PNB,euro$TMI)
plot(euro$URB,euro$TMI)
plot(euro$FEC,euro$TMI)
plot(euro$POP,euro$TMI)Ces figures permettent (avec un peu d’habitude …) d’expliquer les divergences observées entre les coefficients de corrélation de Pearson et de Spearman.
par(mfrow= c(1,2))
X<-euro$PNB
Y<-euro$TMI
corlin<-round(cor(X,Y),3)
reglin<-lm(Y~X)
plot(X,Y, main=paste("Coeff Pearson = ",corlin))
abline(reglin,col="red")
rankX<-rank(X)
rankY<-rank(Y)
corrank<-round(cor(rankX,rankY),3)
regrank<-lm(rankY~rankX)
plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
abline(regrank,col="red")- la relation entre TMI et PNB est clairement non linéaire et forme un arc de parabole plutôt qu’une droite. C’est la raison pour laquelle le coefficient de Spearman était beaucoup plus élevé que celui de Pearson. En effet, lorsque l’on transforme les variables en rang, la courbure disparaît.
par(mfrow= c(1,2))
X<-euro$FEC
Y<-euro$TMI
corlin<-round(cor(X,Y),3)
reglin<-lm(Y~X)
plot(X,Y, main=paste("Coeff Pearson = ",corlin))
abline(reglin,col="red")
rankX<-rank(X)
rankY<-rank(Y)
corrank<-round(cor(rankX,rankY),3)
regrank<-lm(rankY~rankX)
plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
abline(regrank,col="red")- la relation entre TMI et FEC est quant à elle influencée par une valeur exceptionnelle, l’Albanie, qui affiche une fécondité et une mortalité infantile excpetionnnelles. Le point Albanie est situé très en dehors du nuage des autres points et il fait fortement augmenter le coefficient de Pearson. Mais si l’on passe en rang, son effet disparaît et le nuage de point apparaît beaucoup moins linéaire, ce qui explique la valeur plus faible du coefficient de Spearman.
Le point important à retenir (en l’absence de la lecture d’un manuel de statistique plus précis) est qu’on ne doit pas tenter de faire passer une droite dans un nuage de points si celui-ci ne forme pas un ensemble de points régulièrement alignés. Il faut éviter d’utiliser la régression linéaire sur un nuage de points ayant une courbure ou sur un nuage de points comportant des valeurs exceptionnelles. On va voir dans la section suivante quelques “trucs” permettant de corriger la forme d’un nuage de point inadapté à une analyse de régression linéaire (ce qui ne dispense pas de lire un manuel de statistique expliquant ces trucs à l’aide des concepts plus précis d’autocorrélation des résidus, d’hétéroscédasticité, etc…)
2.4.4 Un premier modèle de régression linéaire (quick & dirty …)
Le calcul d’une régression linéaire se fait très facilement avec l’instruction lm qui est l’acronyme de linear model. La particularité de R par rapport à d’autres logiciels est de stocker les résultats d’une régression (et plus généralement d’un modèle statistique quelconque) dans un objet dont on choisit le nom (ici : MonModele). C’est objet est uen sorte de grand sac sans fonds contenant des dizaines de résultats que l’on peut sortir un par un en fonction de ses besoins et surtout de ses connaissances en statistiques… Les débutants se limiteront à demander de sortir du sac à malice un résumé (summary) de la régression qu’ils viennent d’opérer :
MonModele <- lm(euro$TMI~euro$PNB)
summary(MonModele)
Call:
lm(formula = euro$TMI ~ euro$PNB)
Residuals:
Min 1Q Median 3Q Max
-7.921 -3.222 -1.440 1.064 22.344
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.3403638 2.3136443 9.224 5.14e-09 ***
euro$PNB -0.0011403 0.0002647 -4.308 0.000284 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.411 on 22 degrees of freedom
Multiple R-squared: 0.4576, Adjusted R-squared: 0.4329
F-statistic: 18.56 on 1 and 22 DF, p-value: 0.0002843
On trouve dans ce résumé l’essentiel à savoir l’équation de la droite de régression linéaire dans la partie coefficient (TMI = -0.00114 * PNB + 21.34), le pouvoir explicatif du modèle qui est le carré du coefficient de corrélation appelé en anglais R-Squared (0.457 soit 46%) et la significativité générale du modèle qui est fournie par la p-value (0.00028). Cette dernière est en l’occurence conforme au coefficient de corrélation entre TMI et PNB car nous faisons une régression simple. On verra ultérieurement le cas des régressions multiples.
La visualisation de cette droite de régression peut se faire très simplement en traçant le nuage de point avec plot et en ajoutant la droite en utilisant l’instruction abline appliquée au sac à malice où se trouvent les résultats de notre modèle. Même si nous ignorons le contenu détaillé de ce sac, nous savons que R va y trouver les éléments utiles à savoir les coefficients de la droite de régression que nous avons déjà vu dans le résumé. On va faire un petit effort d’habillage en ajoutant à la figure un titre principal (main) et des titres d’axes (xlab,ylab).
<<<<<<< HEADpar(mfrow=c(1,1))
plot(euro$PNB,euro$TMI,
main=" Richesse et mortalité infantile en Europe vers 1988",
ylab="PNB en $ par habitant (1988)",
xlab="Taux de mortalité infantile en p. 1000 (1988)")
abline(MonModele, col="red")par(mfrow=c(1,1))
plot(euro$PNB,euro$TMI,
main=" Richesse et mortalité infantile en Europe vers 1988",
ylab="PNB en $ par habitant (1988)",
xlab="Taux de mortalité infantile en p. 1000 (1988)")
abline(MonModele, col="red")Nous avons réussi à faire notre premier modèle de régression sous R ! Mais la vérité oblige à dire qu’il n’est pas satisfaisant puisque nous avons ajusté une droite dans un nuage de point qui est manifestement de forme non linéaire. Essayons de faire mieux, en montrant que la question statistique peut recouper des question plus théoriques de géographie économique et de géographie politique.
2.4.5 Un deuxième modèle plus satisfaisant sur le plan statistique et théorique
La courbure de la relation entre PNB/habitant et taux de mortalité infantile est certes un problème statistique. Mais elle soulève aussi une question théorique plus intéressante : Est-il raisonnable et logique de penser que plus le PNB/hab. augmente, plus la mortalité infantiel décroît de façon linéaire ?
Notre modèle précédent (Y=aX+b) est en partie absurde puisqu’il nous dit que la mortalité infantile minimale est de 21.34°/°° dans un pays ayant un PNB nul (on sait qu’il existe des mortalité infantile plus forte) et surtout il suppose que chaque fois que le PNB augmente de 1000 $, la mortalité infantile baisse de 1.1 °/°°. Du coup, ce modèle prévoit une absence totale de mortalité infantile pour les pays ayant un revenu supérieur à 15000 $/hab. Et il annonce froidement des mortalités infantile négatives au delà de ce niveau de richesse ! Bref, il y a une inconsistance logique dans un tel modèle, indépendamment du problème statistique de non linéarité.
Or, ce que nous disent les théoriciens du développement comme A. Sen est l’existence de rendements décroissant dans les progrès au fur et à mesure que la richesse augmente. En d’autres termes, la mortalité infantile diminue rapidement quant on passe de l’extrême apuvreté à la pauvreté. Mais ensuite elle baisse de moins en moins vite et son effet devient négligeable au delà d’un seuil de richesse. Le modèle correspondant à ces hypothèses théoriques est donc un modèle de décroissance non linéaire.
Il peut s’agir d’une loi exponentielle négative qui peut s’écrire Y = b.exp(-aX), ce que l’on peut transformer en équation linéaire en le réécrivant sous la forme log(Y) = log(b)+aX. Ou bien d’une loi puissance négative qui peut s’écrire Y = b.X^a et se transforme en équation linéaire log(y) = a.log(X)+log(b). Si l’on n’a pas de préférence théorique pour l’une ou l’autre de ces lois, on peut simplement examiner celle qui donne le meilleur ajustement linéaire après transformation logarithmique :
<<<<<<< HEADpar(mfrow= c(1,2))
X<-euro$PNB
Y<-euro$TMI
corexp<-round(cor(X,log(Y)),3)
regexp<-lm(log(Y)~X)
plot(X,log(Y), main=paste("Modèle exponentiel : r = ",corexp))
abline(regexp,col="red")
corpui<-round(cor(log(X),log(Y)),3)
regpui<-lm(log(Y)~log(X))
plot(log(X),log(Y), main=paste("Modèle puissance : r = ",corpui))
abline(regpui,col="red")par(mfrow= c(1,2))
X<-euro$PNB
Y<-euro$TMI
corexp<-round(cor(X,log(Y)),3)
regexp<-lm(log(Y)~X)
plot(X,log(Y), main=paste("Modèle exponentiel : r = ",corexp))
abline(regexp,col="red")
corpui<-round(cor(log(X),log(Y)),3)
regpui<-lm(log(Y)~log(X))
plot(log(X),log(Y), main=paste("Modèle puissance : r = ",corpui))
abline(regpui,col="red")La comparaison des coefficients de corrélation mais aussi l’examen de la forme des nuages de points plaide clairement en faveur du choix du modèle puissance. Il demeure en effet une nette courbure du nuage de points dans le graphique du modèle exponentiel et l’ajustement obtenu est plus faible. Le modèle puissance affiche un alignement presque parfait des points avec une corrélation très élevée. Donc, si l’on ne dispose pas d’arguments théoriques ou empiriques particuliers en faveur du modèle exponentiel, on retiendra le modèle puissance pour modéliser la relation entre mortalité infantile et PNB/habitant des pays européens en 1988.
2.5 Exercice
Essayez de reproduire la même analyse sur deux autres variables quantitatives de votre choix
3 R-BASE : Variable qualitatives, tableau de contingence, test du chi-2
On suppose que vous avez acquis les premiers réflexes de programmation en R avec les variables quantitatives contenues dans le fichier Europe 1988. VOus y avez appris à expliquer une variable Y quantitative par une variable X pouvant être soit qualitative (analyse de variance, test d’égalité des moyennes), soit quantitative (corrélation, régression)
Dans ce deuxième module, on va essayer d’expliquer une variable qualitative Y par une autre variable qualitative X. On ne traitera pas directement le cas de l’explication d’une variable qualitative Y par une variable quantitative X, car il peut soit se ramener au cas déjà vu dans le cours précédent, soit se résoudre par une transformation de X en une variable X’ qualitative.
L’exemple utilisé ici est fondé sur une véritable enquête de terrain réalisée en 2003 par une étudiante de maitrise de Paris 7, anne-Gaëlle Monnier,à propos de l’électrification d’un village camerounais nommé Yemessoa.
3.1 Définir un environnement de travail
3.1.1 Définir un espace de travail
On vérifie l’emplacement où l’on va trouver les données
# (1.2) Choix du rÈpertoire contenant les donnÈes
monchemin<-"data/yemessoa"
list.files(monchemin)[1] "yemessoa2005.txt" "yemessoa2005.xls"
Nous allons maintenant essayer de charger le fichier de données
3.2 Charger un tableau de données
On utilise ici une procédure adaptée à la lecture de fichiers texte : read.table()
- L’instruction header=TRUE signale que la première ligne donne le nom des variables
- L’instruction *sep = que le séparateur de colonne est une tabulation et non pas un point-virgule
- L’instruction dec=“,” signale que les décimales sont représentées par des virgules et non pas des points comme dans le premier module.
- Linstruction stringsAsFactors = TRUE signifie que les variables de type character seront transformées automatiquement en factor.
3.2.1 Importation d’un tableau
# --------------------------------------------------------
# (2) IMPORTATION ET MISE EN FORME D'UN TABLEAU DE DONNEES
# --------------------------------------------------------
# (2.1) Importation d'un fichier .txt
yem<-read.csv("data/yemessoa/yemessoa2005.txt",
header=TRUE,
sep="\t",
dec=",",
stringsAsFactors = T)On affiche les premières lignes (head) ou les dernières lignes (tail) du tableau pour voir si le chargement s’est bien effectué :
<<<<<<< HEADhead(yem) CODE CON2 CON3 DIST QUAR VIL CLAN SEXE SELF
1 237 Oui 2 3 Elig Messobo Yemessoa II B\xe9ti Homme (1) Petit
2 98 Non 0 20 Nkol Ngu\xe9gu\xe9 Yemessoa I B\xe9loua Homme (1) Petit
3 157 Non 0 30 Biling Bitom Yemessoa I B\xe9loua Homme (1) Petit
4 325 Oui 1 10 Elig Messobo Yemessoa II B\xe9ti Homme (2) Moyen
5 299 Oui 1 10 Nkol Ambenbe Yemessoa III B\xe9ti Homme (3) Grand
BUD TABLE SALON AGE ETU5 ETU2
1 3950 Non Non 19 ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
2 3700 Oui Non 20 ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
3 1250 Non Non 21 ETU5 : BEPC et + ETU2 : Sup\xe9rieur au CEP
4 2200 Oui Non 21 ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
5 4450 Oui Non 21 ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
tail(yem) CODE CON2 CON3 DIST QUAR VIL CLAN SEXE SELF BUD
294 128 Oui 2 5 Biling Bitom Yemessoa I B\xe9ti Homme (3) Grand 3200
295 71 Non 0 35 Nkol Ambenbe Yemessoa III B\xe9ti Homme (3) Grand 2850
296 176 Oui 2 15 Biling Bitom Yemessoa I B\xe9ti Femme (3) Grand 2675
297 178 Oui 2 10 Biling Bitom Yemessoa I B\xe9ti Homme (3) Grand 3750
298 243 Oui 2 25 Mebang Mengoe Yemessoa I B\xe9ti Homme (3) Grand 8900
TABLE SALON AGE ETU5 ETU2
294 Oui Non 81 ETU1: Aucun ETU1: inf\xe9rieur au CEP
295 Oui Non 83 ETU2 : Primaire sans CEP ETU1: inf\xe9rieur au CEP
296 Oui Oui 87 ETU1: Aucun ETU1: inf\xe9rieur au CEP
297 Oui Oui 88 ETU2 : Primaire sans CEP ETU1: inf\xe9rieur au CEP
298 Oui Oui 90 ETU2 : Primaire sans CEP ETU1: inf\xe9rieur au CEP
=======
head(yem)
CODE CON2 CON3 DIST QUAR VIL CLAN SEXE SELF BUD
1 237 Oui 2 3 Elig Messobo Yemessoa II Béti Homme (1) Petit 3950
2 98 Non 0 20 Nkol Nguégué Yemessoa I Béloua Homme (1) Petit 3700
3 157 Non 0 30 Biling Bitom Yemessoa I Béloua Homme (1) Petit 1250
4 325 Oui 1 10 Elig Messobo Yemessoa II Béti Homme (2) Moyen 2200
5 299 Oui 1 10 Nkol Ambenbe Yemessoa III Béti Homme (3) Grand 4450
TABLE SALON AGE ETU5 ETU2
1 Non Non 19 ETU3 : CEP-CAP ETU2 : Supérieur au CEP
2 Oui Non 20 ETU3 : CEP-CAP ETU2 : Supérieur au CEP
3 Non Non 21 ETU5 : BEPC et + ETU2 : Supérieur au CEP
4 Oui Non 21 ETU3 : CEP-CAP ETU2 : Supérieur au CEP
5 Oui Non 21 ETU3 : CEP-CAP ETU2 : Supérieur au CEP
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
tail(yem)
CODE CON2 CON3 DIST QUAR VIL CLAN SEXE SELF BUD
294 128 Oui 2 5 Biling Bitom Yemessoa I Béti Homme (3) Grand 3200
295 71 Non 0 35 Nkol Ambenbe Yemessoa III Béti Homme (3) Grand 2850
296 176 Oui 2 15 Biling Bitom Yemessoa I Béti Femme (3) Grand 2675
297 178 Oui 2 10 Biling Bitom Yemessoa I Béti Homme (3) Grand 3750
298 243 Oui 2 25 Mebang Mengoe Yemessoa I Béti Homme (3) Grand 8900
TABLE SALON AGE ETU5 ETU2
294 Oui Non 81 ETU1: Aucun ETU1: inférieur au CEP
295 Oui Non 83 ETU2 : Primaire sans CEP ETU1: inférieur au CEP
296 Oui Oui 87 ETU1: Aucun ETU1: inférieur au CEP
297 Oui Oui 88 ETU2 : Primaire sans CEP ETU1: inférieur au CEP
298 Oui Oui 90 ETU2 : Primaire sans CEP ETU1: inférieur au CEP
>>>>>>> d676e259c8be1e9f88c818dd0e1c0cc0a591abdd
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
3.2.2 Vérification du type des variables
On a vu qu’on peut s’éviter beaucoup d’ennuis en contrôlant le type des variables qui a été attribué aux différentes colonnes par R lors de la lecture d’un fichier. Pour cela on commence par regarder quels sont les types de variables du tableau à l’aide de l’instruction str() :
str(yem)'data.frame': 299 obs. of 15 variables:
$ CODE : int 237 98 157 325 299 91 215 48 50 154 ...
$ CON2 : Factor w/ 2 levels "Non","Oui": 2 1 1 2 2 1 2 2 2 1 ...
$ CON3 : int 2 0 0 1 1 0 1 2 2 0 ...
$ DIST : int 3 20 30 10 10 1000 20 15 30 50 ...
$ QUAR : Factor w/ 9 levels "Biling Bitom",..: 4 8 1 4 6 7 1 4 4 1 ...
$ VIL : Factor w/ 3 levels "Yemessoa I","Yemessoa II",..: 2 1 1 2 3 3 1 2 2 1 ...
$ CLAN : Factor w/ 2 levels "B\xe9loua","B\xe9ti": 2 1 1 2 2 2 2 2 2 2 ...
$ SEXE : Factor w/ 2 levels "Femme","Homme": 2 2 2 2 2 1 2 2 2 2 ...
$ SELF : Factor w/ 3 levels "(1) Petit","(2) Moyen",..: 1 1 1 2 3 1 2 2 3 1 ...
$ BUD : int 3950 3700 1250 2200 4450 1500 1150 8550 13575 1750 ...
$ TABLE: Factor w/ 2 levels "Non","Oui": 1 2 1 2 2 2 2 1 1 1 ...
$ SALON: Factor w/ 2 levels "Non","Oui": 1 1 1 1 1 2 1 1 1 1 ...
$ AGE : int 19 20 21 21 21 22 24 25 26 28 ...
$ ETU5 : Factor w/ 5 levels "ETU1: Aucun",..: 3 3 5 3 3 2 5 2 5 2 ...
$ ETU2 : Factor w/ 2 levels "ETU1: inf\xe9rieur au CEP",..: 2 2 2 2 2 1 2 1 2 1 ...
Il y a visiblement des erreurs de codage que l’on va corriger.
yem$CODE<-as.character(yem$CODE)
yem$CON3<-as.factor(yem$CON3)
levels (yem$CON3)[1] "0" "1" "2"
levels(yem$CON3)<-c("(0) : non connecté","(1) informel","(2) officiel")
yem$DIST<-as.numeric(yem$DIST)
str(yem)'data.frame': 299 obs. of 15 variables:
$ CODE : chr "237" "98" "157" "325" ...
$ CON2 : Factor w/ 2 levels "Non","Oui": 2 1 1 2 2 1 2 2 2 1 ...
$ CON3 : Factor w/ 3 levels "(0) : non connecté",..: 3 1 1 2 2 1 2 3 3 1 ...
$ DIST : num 3 20 30 10 10 1000 20 15 30 50 ...
$ QUAR : Factor w/ 9 levels "Biling Bitom",..: 4 8 1 4 6 7 1 4 4 1 ...
$ VIL : Factor w/ 3 levels "Yemessoa I","Yemessoa II",..: 2 1 1 2 3 3 1 2 2 1 ...
$ CLAN : Factor w/ 2 levels "B\xe9loua","B\xe9ti": 2 1 1 2 2 2 2 2 2 2 ...
$ SEXE : Factor w/ 2 levels "Femme","Homme": 2 2 2 2 2 1 2 2 2 2 ...
$ SELF : Factor w/ 3 levels "(1) Petit","(2) Moyen",..: 1 1 1 2 3 1 2 2 3 1 ...
$ BUD : int 3950 3700 1250 2200 4450 1500 1150 8550 13575 1750 ...
$ TABLE: Factor w/ 2 levels "Non","Oui": 1 2 1 2 2 2 2 1 1 1 ...
$ SALON: Factor w/ 2 levels "Non","Oui": 1 1 1 1 1 2 1 1 1 1 ...
$ AGE : int 19 20 21 21 21 22 24 25 26 28 ...
$ ETU5 : Factor w/ 5 levels "ETU1: Aucun",..: 3 3 5 3 3 2 5 2 5 2 ...
$ ETU2 : Factor w/ 2 levels "ETU1: inf\xe9rieur au CEP",..: 2 2 2 2 2 1 2 1 2 1 ...
Et voilà : désormais notre tableau est désormais cohérent. On pourra encore l’améliorer par la suite en transformant des variables quantitatives en variables qualitatives
3.2.3 Un petit résumé statistique
La récompense d’une bonne définition du tableau est d’obtenir pour chaque variable un résumé adapté à son type
summary(yem) CODE CON2 CON3 DIST
Length:299 Non:165 (0) : non connecté:165 Min. : 1.0
Class :character Oui:134 (1) informel : 49 1st Qu.: 15.0
Mode :character (2) officiel : 85 Median : 25.0
Mean : 370.4
3rd Qu.: 500.0
QUAR VIL CLAN SEXE
Biling Bitom :61 Yemessoa I :150 B\xe9loua: 32 Femme: 51
Nkol Ambenbe :60 Yemessoa II : 58 B\xe9ti :267 Homme:248
Elig Messobo :57 Yemessoa III: 91
Mebang Mengoe :30
Nkol Ngu\xe9gu\xe9:30
SELF BUD TABLE SALON AGE
(1) Petit: 54 Min. : 350 Non: 86 Non:191 Min. :19.00
(2) Moyen: 90 1st Qu.: 2200 Oui:213 Oui:108 1st Qu.:40.00
(3) Grand:155 Median : 3300 Median :50.00
Mean : 4937 Mean :51.14
3rd Qu.: 5700 3rd Qu.:63.00
ETU5 ETU2
ETU1: Aucun : 52 ETU1: inf\xe9rieur au CEP :164
ETU2 : Primaire sans CEP:112 ETU2 : Sup\xe9rieur au CEP:135
ETU3 : CEP-CAP : 72
ETU4 : Coll\xe8ge : 34
ETU5 : BEPC et + : 29
[ getOption("max.print") est atteint -- 2 lignes omises ]
Avant d’aller plus loin, commentez rapidement ce tableau pour bien prendre connaissance de chacune des variables.
3.3 Mise en relation de deux variables qualitatives
Le niveau d’électrification des ménages (Y) est-il lié au genre du chef de ménage (X) ?
Notre objectif est le même que dans le module précédent : chercher s’il existe une relation significative entre X et Y puis en décrire la forme.
3.3.1 Analyse de la variable dépendante (Y)
Deux petites commandes pour dénombrer l’effectif des classes et les représenter
# Analyse de la variable dépendante
Y<-yem$CON2
#calcul des fréquences
summary(Y)Non Oui
165 134
<<<<<<< HEAD
#diagramme en bâtons
plot(Y) Commentaire : …
#diagramme en bâtons
plot(Y) Commentaire : …
3.3.2 Analyse de la variable indépendante (X)
Même programme …
# Analyse de la variable dépendante
X<-yem$SEXE
#calcul des fréquences
summary(X)Femme Homme
51 248
<<<<<<< HEAD
#diagramme en bâtons
plot(X)Commentaire : …
2.3.3 Visualisation de la distribution croisée de Y et X
plot(Y~X)#diagramme en bâtons
plot(X)Commentaire : …
3.3.3 Visualisation de la distribution croisée de Y et X
plot(Y~X)Commentaire : …
3.3.4 Calcul du tableau de contingence et des trois tableaux de pourcentage associés
tabcont<-table(Y,X)
tabcont X
Y Femme Homme
Non 36 129
Oui 15 119
prop.table(tabcont,margin=) X
Y Femme Homme
Non 0.12040134 0.43143813
Oui 0.05016722 0.39799331
prop.table(tabcont,margin=1) X
Y Femme Homme
Non 0.2181818 0.7818182
Oui 0.1119403 0.8880597
prop.table(tabcont,margin=2) X
Y Femme Homme
Non 0.7058824 0.5201613
Oui 0.2941176 0.4798387
Commentaire : …
3.3.5 Calcul du tableau théorique et test du chi-2
test<-chisq.test(tabcont)
test$observed X
Y Femme Homme
Non 36 129
Oui 15 119
test$expected X
Y Femme Homme
Non 28.14381 136.8562
Oui 22.85619 111.1438
test$residuals X
Y Femme Homme
Non 1.4808817 -0.6715519
Oui -1.6432738 0.7451937
test
Pearson's Chi-squared test with Yates' continuity correction
data: tabcont
X-squared = 5.1726, df = 1, p-value = 0.02295
Commentaire : ***
2.3.6 Une jolie figure pour interpréter la relation (si elle est significative)
mosaicplot(~Y+X,
main = "Relation entre X et Y", shade = TRUE, legend = TRUE)3.3.6 Une jolie figure pour interpréter la relation (si elle est significative)
mosaicplot(~Y+X,
main = "Relation entre X et Y", shade = TRUE, legend = TRUE)3.4 Exercice
Essayez de reproduire la même analyse sur deux autres variables qualitatives de votre choix.
Bibliographie
Annexes
Infos session
| setting | value |
|---|---|
| version | R version 4.0.2 (2020-06-22) |
| os | macOS Catalina 10.15.7 |
| system | x86_64, darwin17.0 |
| ui | X11 |
| language | (EN) |
| collate | fr_FR.UTF-8 |
| ctype | fr_FR.UTF-8 |
| tz | Europe/Paris |
| date | 2021-09-22 |
| package | ondiskversion | source |
|---|---|---|
| knitr | 1.34 | CRAN (R 4.0.2) |
| rmarkdown | 2.11 | CRAN (R 4.0.2) |
| rzine | 0.1.0 | gitlab (rzine/package@a94bf55) |
Citation
@Manual{ficheRzine,
title = {Titre de la fiche},
author = {{Auteur.e.s}},
organization = {Rzine},
year = {202x},
url = {http://rzine.fr/},
}